
* -------------------
* Title: 3_analysis.do
* -------------------

* Preface: This is the THIRD of four code files to reproduce all results
* reported in Joo, Elwert, and Munk 2024 "Labor Market Consequences of 
* Grandparenthood" from source data stored on Statistics Denmark servers.
* See ReadMe for details on data access.  This do file has been minimally
* redacted by Statistics Denmark staff to remove identifying information
* in compliance with applicable law. 

* Content: this code executes all statistical models reported in the paper.



* Set up
* ---------------
version 17.0
set seed 12345
set more off
set type float
set matsize 11000

* Define macros
* ---------------
global home "E:\workdata\704121\Wontak"
global raw "E:\workdata\704121\Wontak\raw"
global data "E:\workdata\704121\Wontak\data"
global result "E:\workdata\704121\Wontak\result"

global y "g1_empd g1_empdi g1_incomel"
global d1 "d"
global dc1 "dc"
global f1 "female teenb_g1 ager1 ager2 year yearr t1 g"
global g1 "births_t g1_workexp g1_empd g1_empdi g1_incomel g1_yed g1_ced g1_psych g1_ccs g1_dhosp g1_mstatr g1_partner p_empdi p_incomeli p_psych p_ccs p_dhosp"
global g21 "g2_1_partner g2_1_workexp g2_1_empd g2_1_incomel g2_1_yed g2_1_ced g2_1_psych g2_1_ccs g2_1_dhosp"
global g2 "g2 g3 g2_female g2_partner g2_empd g2_incomel g2_yed g2_ced g2_psych g2_ccs g2_dhosp"
global h "h_incomeg hhsize hh_g2 hh_u18 h_welfare"

global v1_g1 "L1.c.births_t L1.c.g1_workexp L1.c.g1_empd L1.i.g1_empdi L1.c.g1_incomel L1.c.g1_yed L1.i.g1_ced L1.c.g1_psych L1.c.g1_ccs L1.c.g1_dhosp L1.i.g1_partner L1.i.p_empdi L1.i.p_incomeli L1.i.p_psych L1.i.p_ccs L1.i.p_dhosp"
global v2_g1 "L2.c.births_t L2.c.g1_workexp L2.c.g1_empd L2.i.g1_empdi L2.c.g1_incomel L2.c.g1_yed L2.i.g1_ced L2.c.g1_psych L2.c.g1_ccs L2.c.g1_dhosp L2.i.g1_partner L2.i.p_empdi L2.i.p_incomeli L2.i.p_psych L2.i.p_ccs L2.i.p_dhosp"
global v1_g21 "L1.i.g2_1_partner L1.c.g2_1_workexp L1.c.g2_1_empd L1.c.g2_1_incomel L1.c.g2_1_yed L1.i.g2_1_ced L1.c.g2_1_psych L1.c.g2_1_ccs L1.c.g2_1_dhosp"
global v2_g21 "L2.i.g2_1_partner L2.c.g2_1_workexp L2.c.g2_1_empd L2.c.g2_1_incomel L2.c.g2_1_yed L2.i.g2_1_ced L2.c.g2_1_psych L2.c.g2_1_ccs L2.c.g2_1_dhosp"
global v1_g2 "L1.c.g2 L1.c.g3 L1.c.g2_female L1.c.g2_partner L1.c.g2_empd L1.c.g2_incomel L1.c.g2_yed L1.c.g2_ced L1.c.g2_psych L1.c.g2_ccs L1.c.g2_dhosp"
global v2_g2 "L2.c.g2 L2.c.g3 L2.c.g2_female L2.c.g2_partner L2.c.g2_empd L2.c.g2_incomel L2.c.g2_yed L2.c.g2_ced L2.c.g2_psych L2.c.g2_ccs L2.c.g2_dhosp"
global v1_h "L1.c.h_incomeg L1.c.hhsize L1.c.hh_g2 L1.c.hh_u18 L1.c.h_welfare"
global v2_h "L2.c.h_incomeg L2.c.hhsize L2.c.hh_g2 L2.c.hh_u18 L2.c.h_welfare"

global t1_g1 "c.births_t_t1 c.g1_workexp_t1 c.g1_empd_t1 i.g1_empdi_t1 c.g1_incomel_t1 c.g1_yed_t1 c.g1_psych_t1 c.g1_ccs_t1 c.g1_dhosp_t1 i.g1_mstatr_t1 i.p_empdi_t1 i.p_incomeli_t1 i.p_psych_t1 i.p_ccs_t1 i.p_dhosp_t1"
global t1_g21 "c.g2_1_yed_t1 c.g2_1_psych_t1 c.g2_1_ccs_t1 c.g2_1_dhosp_t1"
global t1_g2 "c.g2_t1 c.g2_female_t1"
global t1_h "c.h_incomeg_t1 c.hhsize_t1 c.hh_g2_t1 c.hh_u18_t1 c.h_welfare_t1"

global f1a "i.ager1 i.year i.t1"
global f1b "i.female i.teenb_g1 i.g"

global v1_g1e "L1.c.births_t L1.c.g1_workexp L1.c.g1_incomel L1.c.g1_yed L1.i.g1_ced L1.c.g1_psych L1.c.g1_ccs L1.c.g1_dhosp L1.i.g1_partner L1.i.p_empdi L1.i.p_incomeli L1.i.p_psych L1.i.p_ccs L1.i.p_dhosp"
global v1_g1i "L1.c.births_t L1.c.g1_workexp L1.c.g1_empd L1.i.g1_empdi L1.c.g1_yed L1.i.g1_ced L1.c.g1_psych L1.c.g1_ccs L1.c.g1_dhosp L1.i.g1_partner L1.i.p_empdi L1.i.p_incomeli L1.i.p_psych L1.i.p_ccs L1.i.p_dhosp"

cd $home


* Load and prep data
* -----------------------

use $data\msm_230115.dta, clear

* dcr dcr_0 dcr_1
replace g1_empdi=100*g1_empdi

gen d_0=0 if s==1
replace d_0=d if tb==0
gen d_1=0 if s==1
replace d_1=d if tb==1

recode dc (11/max=11), gen(dcr)
gen dcr_0=0 if s==1
replace dcr_0=dcr if tb==0
gen dcr_1=0 if s==1
replace dcr_1=dcr if tb==1

* define subsamples: rich (individual labor income) / hrich (household income)
_pctile g1_incomel_b3 if t1==3 & female==0, p(33.33333, 66.66667)
disp r(r1)
disp r(r2)
local t1=r(r1)
local t2=r(r2)
recode g1_incomel_b3 (min/`t1'=1)(`t1'/`t2'=2)(`t2'/max=3) if female==0, gen(richm)
_pctile g1_incomel_b3 if t1==3 & female==1, p(33.33333, 66.66667)
disp r(r1)
disp r(r2)
local t1=r(r1)
local t2=r(r2)
recode g1_incomel_b3 (min/`t1'=1)(`t1'/`t2'=2)(`t2'/max=3) if female==1, gen(richf)
gen rich=richm
replace rich=richf if female==1
drop richm richf

_pctile h_incomeg_b3 if t1==3 & female==0, p(33.33333, 66.66667)
disp r(r1)
disp r(r2)
local t1=r(r1)
local t2=r(r2)
recode h_incomeg_b3 (min/`t1'=1)(`t1'/`t2'=2)(`t2'/max=3) if female==0, gen(hrichm)
_pctile h_incomeg_b3 if t1==3 & female==1, p(33.33333, 66.66667)
disp r(r1)
disp r(r2)
local t1=r(r1)
local t2=r(r2)
recode h_incomeg_b3 (min/`t1'=1)(`t1'/`t2'=2)(`t2'/max=3) if female==1, gen(hrichf)
gen hrich=hrichm
replace hrich=hrichf if female==1
drop hrichm hrichf

log close

* matrices
foreach i in empd empdi incomel{
	forvalues j=0/1{
		mat `i'`j'_b=J(1,7,.)
		mat `i'`j'=J(11,7,.)
		foreach s in rich hrich{
			forvalues a=1/3{
				mat `i'`j'_`s'`a'=J(11,7,.)
			}
		}
		forvalues k=0/1{
			mat `i'`j'`k'_b=J(1,7,.)
			mat `i'`j'`k'=J(11,7,.)
			foreach s in rich hrich{
				forvalues a=1/3{
					mat `i'`j'`k'_`s'`a'=J(11,7,.)
				}
			}
		}
	}
}

********************************************************************************
***************************************************** Inverse probability weight 
********************************************************************************

* Baseline model (binary) (Table 1, A4, A5, and A6)
* -----------------------------------------------------

foreach i in empd empdi incomel{
	forvalues j=0/1{
		eststo: reg g1_`i' $f1a i.d i.tb#i.d if t1>=3 & s==1 & female==`j', vce(cl pnrr)
		eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.d i.tb#i.d if t1>=3 & s==1 & female==`j', vce(cl pnrr)
		eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h $v1_g1 $v1_g21 $v1_g2 $v1_h) i.d i.tb#i.d if t1>=3 & s==1 & female==`j', vce(cl pnrr)
		if `i'=="empd"|`i'=="empdi"{
			eststo: areg g1_`i' $f1a $v1_g1e $v1_g21 $v1_g2 $v1_h i.$d1 i.tb#i.$d1 if t1>=3 & s==1 & female==`j', absorb(pnrr) vce(cl pnrr)
		}
		if `i'=="incomel"{
			eststo: areg g1_`i' $f1a $v1_g1i $v1_g21 $v1_g2 $v1_h i.$d1 i.tb#i.$d1 if t1>=3 & s==1 & female==`j', absorb(pnrr) vce(cl pnrr)
		}
		eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.d i.tb#i.d [pw=w] if t1>=3 & s==1 & female==`j', vce(cl pnrr)
		mat `i'`j'_b[1,7]=2*ttail(e(df_r),abs(_b[1.tb#1.d]/_se[1.tb#1.d]))
	}
	forvalues j=0/1{
		forvalues k=0/1{
			eststo: reg g1_`i' $f1a i.d i.tb#i.d if t1>=3 & s==1 & female==`j' & g==`k', vce(cl pnrr)
			eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.d i.tb#i.d if t1>=3 & s==1 & female==`j' & g==`k', vce(cl pnrr)
			eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h $v1_g1 $v1_g21 $v1_g2 $v1_h) i.d i.tb#i.d if t1>=3 & s==1 & female==`j' & g==`k', vce(cl pnrr)
			if `i'=="empd"|`i'=="empdi"{
				eststo: areg g1_`i' $f1a $v1_g1e $v1_g21 $v1_g2 $v1_h i.$d1 i.tb#i.$d1 if t1>=3 & s==1 & female==`j' & g==`k', absorb(pnrr) vce(cl pnrr)
			}
			if `i'=="incomel"{
				eststo: areg g1_`i' $f1a $v1_g1i $v1_g21 $v1_g2 $v1_h i.$d1 i.tb#i.$d1 if t1>=3 & s==1 & female==`j' & g==`k', absorb(pnrr) vce(cl pnrr)
			}
			eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.d i.tb#i.d [pw=w] if t1>=3 & s==1 & female==`j' & g==`k', vce(cl pnrr)
			mat `i'`j'`k'_b[1,7]=2*ttail(e(df_r),abs(_b[1.tb#1.d]/_se[1.tb#1.d]))
		}
	}
	estout * using $result\msm_`i'_binary.txt, replace noomit nobase cells(b(star fmt(3)) se(par fmt(3)) p(par fmt(3))) stats(N F r2)
	estimate clear
}

foreach i in empd empdi incomel{
	forvalues j=0/1{
		reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.d_0 i.d_1 [pw=w] if t1>=3 & s==1 & female==`j', vce(cl pnrr)
		capture mat `i'`j'_b[1,1]=_b[1.d_0]
		capture mat `i'`j'_b[1,2]=_b[1.d_0]-_se[1.d_0]*invttail(e(df_r),0.005)
		capture mat `i'`j'_b[1,3]=_b[1.d_0]+_se[1.d_0]*invttail(e(df_r),0.005)
		capture mat `i'`j'_b[1,4]=_b[1.d_1]
		capture mat `i'`j'_b[1,5]=_b[1.d_1]-_se[1.d_1]*invttail(e(df_r),0.005)
		capture mat `i'`j'_b[1,6]=_b[1.d_1]+_se[1.d_1]*invttail(e(df_r),0.005)
	}
	forvalues j=0/1{
		forvalues k=0/1{
			reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.d_0 i.d_1 [pw=w] if t1>=3 & s==1 & female==`j' & g==`k', vce(cl pnrr)
			capture mat `i'`j'`k'_b[1,1]=_b[1.d_0]
			capture mat `i'`j'`k'_b[1,2]=_b[1.d_0]-_se[1.d_0]*invttail(e(df_r),0.005)
			capture mat `i'`j'`k'_b[1,3]=_b[1.d_0]+_se[1.d_0]*invttail(e(df_r),0.005)
			capture mat `i'`j'`k'_b[1,4]=_b[1.d_1]
			capture mat `i'`j'`k'_b[1,5]=_b[1.d_1]-_se[1.d_1]*invttail(e(df_r),0.005)
			capture mat `i'`j'`k'_b[1,6]=_b[1.d_1]+_se[1.d_1]*invttail(e(df_r),0.005)
		}
	}
}

* Trajectory of grandparenthood effects
* ---------------------------------------

foreach i in empd empdi incomel{
	forvalues j=0/1{
		eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.dcr i.tb#i.dcr [pw=w] if t1>=3 & s==1 & female==`j', vce(cl pnrr)
		forvalues r=1/11{
			mat `i'`j'[`r',7]=2*ttail(e(df_r),abs(_b[1.tb#`r'.dcr]/_se[1.tb#`r'.dcr]))
		}
	}
	forvalues j=0/1{
		forvalues k=0/1{
			eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.dcr i.tb#i.dcr [pw=w] if t1>=3 & s==1 & female==`j' & g==`k', vce(cl pnrr)
			forvalues r=1/11{
				mat `i'`j'`k'[`r',7]=2*ttail(e(df_r),abs(_b[1.tb#`r'.dcr]/_se[1.tb#`r'.dcr]))
			}
		}
	}
	estout * using $result\msm_`i'_nonlinear.txt, replace noomit nobase cells(b(star fmt(3)) se(par fmt(3)) p(par fmt(3))) stats(N F r2)
	estimate clear
}

foreach i in empd empdi incomel{
	forvalues j=0/1{
		reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.dcr_0 i.dcr_1 [pw=w] if t1>=3 & s==1 & female==`j', vce(cl pnrr)
		forvalues r=1/11{
			capture mat `i'`j'[`r',1]=_b[`r'.dcr_0]
			capture mat `i'`j'[`r',2]=_b[`r'.dcr_0]-_se[`r'.dcr_0]*invttail(e(df_r),0.005)
			capture mat `i'`j'[`r',3]=_b[`r'.dcr_0]+_se[`r'.dcr_0]*invttail(e(df_r),0.005)
			capture mat `i'`j'[`r',4]=_b[`r'.dcr_1]
			capture mat `i'`j'[`r',5]=_b[`r'.dcr_1]-_se[`r'.dcr_1]*invttail(e(df_r),0.005)
			capture mat `i'`j'[`r',6]=_b[`r'.dcr_1]+_se[`r'.dcr_1]*invttail(e(df_r),0.005)
		}
	}
	forvalues j=0/1{
		forvalues k=0/1{
			reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.dcr_0 i.dcr_1 [pw=w] if t1>=3 & s==1 & female==`j' & g==`k', vce(cl pnrr)
			forvalues r=1/11{
				capture mat `i'`j'`k'[`r',1]=_b[`r'.dcr_0]
				capture mat `i'`j'`k'[`r',2]=_b[`r'.dcr_0]-_se[`r'.dcr_0]*invttail(e(df_r),0.005)
				capture mat `i'`j'`k'[`r',3]=_b[`r'.dcr_0]+_se[`r'.dcr_0]*invttail(e(df_r),0.005)
				capture mat `i'`j'`k'[`r',4]=_b[`r'.dcr_1]
				capture mat `i'`j'`k'[`r',5]=_b[`r'.dcr_1]-_se[`r'.dcr_1]*invttail(e(df_r),0.005)
				capture mat `i'`j'`k'[`r',6]=_b[`r'.dcr_1]+_se[`r'.dcr_1]*invttail(e(df_r),0.005)
			}
		}
	}
}

* Trajectory of grandparenthood effects by individual labor income / household income
* -------------------------------------------------------------------------------------

foreach i in empd empdi incomel{
	foreach s in rich hrich{
		forvalues a=1/3{
			forvalues j=0/1{
				eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.dcr i.tb#i.dcr [pw=w] if t1>=3 & s==1 & female==`j' & `s'==`a', vce(cl pnrr)
				forvalues r=1/11{
					mat `i'`j'_`s'`a'[`r',7]=2*ttail(e(df_r),abs(_b[1.tb#`r'.dcr]/_se[1.tb#`r'.dcr]))
				}
			}
			forvalues j=0/1{
				forvalues k=0/1{
					eststo: reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.dcr i.tb#i.dcr [pw=w] if t1>=3 & s==1 & female==`j' & g==`k' & `s'==`a', vce(cl pnrr)
					forvalues r=1/11{
						mat `i'`j'`k'_`s'`a'[`r',7]=2*ttail(e(df_r),abs(_b[1.tb#`r'.dcr]/_se[1.tb#`r'.dcr]))
					}
				}
			}
		}
		estout * using $result\msm_`i'_nonlinear_`s'.txt, replace noomit nobase cells(b(star fmt(3)) se(par fmt(3)) p(par fmt(3))) stats(N F r2)
		estimate clear
	}
}

foreach i in empd empdi incomel{
	foreach s in rich hrich{
		forvalues a=1/3{
			forvalues j=0/1{
				reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.dcr_0 i.dcr_1 [pw=w] if t1>=3 & s==1 & female==`j' & `s'==`a', vce(cl pnrr)
				forvalues r=1/11{
					capture mat `i'`j'_`s'`a'[`r',1]=_b[`r'.dcr_0]
					capture mat `i'`j'_`s'`a'[`r',2]=_b[`r'.dcr_0]-_se[`r'.dcr_0]*invttail(e(df_r),0.005)
					capture mat `i'`j'_`s'`a'[`r',3]=_b[`r'.dcr_0]+_se[`r'.dcr_0]*invttail(e(df_r),0.005)
					capture mat `i'`j'_`s'`a'[`r',4]=_b[`r'.dcr_1]
					capture mat `i'`j'_`s'`a'[`r',5]=_b[`r'.dcr_1]-_se[`r'.dcr_1]*invttail(e(df_r),0.005)
					capture mat `i'`j'_`s'`a'[`r',6]=_b[`r'.dcr_1]+_se[`r'.dcr_1]*invttail(e(df_r),0.005)
				}
			}
			forvalues j=0/1{
				forvalues k=0/1{
					reg g1_`i' $f1a (c.t1##c.t1)##($f1b $t1_g1 $t1_g21 $t1_g2 $t1_h) i.dcr_0 i.dcr_1 [pw=w] if t1>=3 & s==1 & female==`j' & g==`k' & `s'==`a', vce(cl pnrr)
					forvalues r=1/11{
						capture mat `i'`j'`k'_`s'`a'[`r',1]=_b[`r'.dcr_0]
						capture mat `i'`j'`k'_`s'`a'[`r',2]=_b[`r'.dcr_0]-_se[`r'.dcr_0]*invttail(e(df_r),0.005)
						capture mat `i'`j'`k'_`s'`a'[`r',3]=_b[`r'.dcr_0]+_se[`r'.dcr_0]*invttail(e(df_r),0.005)
						capture mat `i'`j'`k'_`s'`a'[`r',4]=_b[`r'.dcr_1]
						capture mat `i'`j'`k'_`s'`a'[`r',5]=_b[`r'.dcr_1]-_se[`r'.dcr_1]*invttail(e(df_r),0.005)
						capture mat `i'`j'`k'_`s'`a'[`r',6]=_b[`r'.dcr_1]+_se[`r'.dcr_1]*invttail(e(df_r),0.005)
					}
				}
			}
		}
	}
}

* Export matrices
* ----------------*

clear
foreach i in empd empdi incomel{
	forvalues j=0/1{
		svmat `i'`j'_b, n(`i'`j'_b)
		forvalues k=0/1{
			svmat `i'`j'`k'_b, n(`i'`j'`k'_b)
		}
	}
	forvalues j=0/1{
		svmat `i'`j', n(`i'`j'_o)
		forvalues k=0/1{
			svmat `i'`j'`k', n(`i'`j'`k'_o)
		}
	}
	foreach s in rich hrich{
		forvalues a=1/3{
			forvalues j=0/1{
				svmat `i'`j'_`s'`a', n(`i'`j'_`s'`a'_o)
				forvalues k=0/1{
					svmat `i'`j'`k'_`s'`a', n(`i'`j'`k'_`s'`a'_o)
				}
			}
		}
	}
}
save "$result\graph_msm.dta", replace
